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STATISTICAL CHALLENGES IN THE ANALYSIS OF COSMIC 
MICROWAVE BACKGROUND RADIATION 
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C^ ■ An enormous amount of observations on Cosmic Microwave Back- 

ground radiation has been collected in the last decade, and much 
more data are expected in the near future from planned or operating 
l/") ' satellite missions. These datasets are a goldmine of information for 

Cosmology and Theoretical Physics; their efficient exploitation posits 
several intriguing challenges from the statistical point of view. In this 
paper we review a number of open problems in CMB data analysis 
and we present applications to observations from the WMAP mission. 
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1. Introduction. 



1.1. Cosmological background. Cosmology is now developing into a ma- 
*>. , ture observational science, with a vast array of different experiments that 

\£) \ yield datasets of astonishing magnitude and nearly as great challenges for 

theoretical and applied statisticians. Datasets are now available on a large 
variety of different phenomena, but the leading part in cosmological research 
has been played over the last 15 years by the analysis of Cosmic Microwave 
Background (CMB) radiation, an area which has already led to Nobel Prizes 
for Physics in 1978 and in 2006. 

The nature of CMB can be loosely explained as follows [see, e.g., Dodelson 
(2003) for a textbook account]. According to the standard cosmological 
model, the Universe that we currently observe originated approximately 



H \ 13.7 billion years ago in a very hot and dense state, in what of course is 

universally known as the Big Bang. Neglecting fundamental physics in the 
first fractions of seconds, we can naively imagine a fluid state where matter 
was completely ionized, that is, the kinetic energy of electrons was much 
stronger than the electrical attraction of protons, so that no stable atomic 
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nuclei could form. It is a consequence of quantum principles that a free elec- 
tron has a much larger cross-section than when it is bound in a nucleus; 
loosely speaking, as a consequence, the probability of interactions between 
photons and electrons is so high that the mean free path of the former was 
very short and the Universe was consequently "opaque." As the Universe 
expands, the mean energy content decreases, that is, the fluid of matter and 
radiation cools down; the mean kinetic energy of the electrons decreases 
as well until it reaches a critical value where it is no longer sufficient to 
compensate the electromagnetic attraction of the protons; stable (and neu- 
tral) hydrogen atoms are then formed. This change of state occurs at the 
so-called "age of recombination," which is currently reckoned to have taken 
place 3.7 x 10 5 years after the Big Bang, that is, when the Universe had only 
the 0.003% of its current age. At the age of recombination, the probability 
of interactions became so small that, as a first approximation, photons could 
start to travel freely. Neglecting second order effects, we can assume they 
had no further interaction up to the present epoch. 

The remarkable consequence of this mechanism is that the Universe is 
embedded in a uniform radiation that provides pictures of its state nearly 
1.37 x 10 10 years ago; this is exactly the above-mentioned CMB radiation. 
The existence of CMB was predicted by G.Gamow in a series of papers 
in the forties; it was later discovered fortuitously by Penzias and Wilson in 
1965 — for this discovery they earned the Nobel Prize for Physics in 1978. For 
several years further experiments were only able to confirm the existence of 
the radiation, and to test its adherence to the Planckian curve of blackbody 
emission, as predicted by theorists. A major breakthrough occurred with 
NASA satellite mission COBE, which was launched in 1989 and publicly 
released the first full-sky maps of radiation in 1992; for these maps Smoot 
and Mather earned the Noble Prize for Physics in 2006 [Smoot et al. (1992)]. 

The nature of these maps deserves further explanation. CMB is dis- 
tributed in remarkably uniform fashion over the sky, with deviations in the 
order of 10 -4 with respect to the mean value (corresponding to 2.731 Kelvin 
degrees). The attempts to understand this uniformity have led to very impor- 
tant developments in cosmology, primarily the inflationary scenario which 
now dominates the theoretical landscape. Even more important, though, are 
the tiny fluctuations around this mean value, which provided the seeds for 
stars and galaxies to form out of gravitational instability. Measuring and 
understanding the nature of these fluctuations has then been the core of an 
enormous amount of experimental and theoretical research. In particular, 
their stochastic properties yield a goldmine of information on a variety of 
extremely important issues on astrophysics and cosmology, and on many 
problems at the frontier of fundamental physics. 

To mention just a few of these problems, we recall the issues concerning 
the matter content of the Universe, its global geometry, the existence and 
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nature of (nonbaryonic) dark matter, the existence and nature of dark en- 
ergy, which is related to Einstein's cosmological constant, and many others. 
The next experimental landmark in CMB analysis followed in 2000, when 
two balloon-borne experiments, BOOMERANG and MAXIMA, yielded the 
first high-resolution observations on small patches of the sky (less than 10° 
squared). These observations led to the first constraints on the global ge- 
ometry of the Universe, which was found to be (very close to) Euclidean. 
Another major breakthrough followed with the 2003, 2007 and 2008 data 
releases from another NASA satellite experiment, WMAP (the data are 
publicly available on the web site http://lambda.gsfc.nasa.gov/). Such 
data releases yielded measurement of the correlation structure of the random 
field up to a resolution of about 0.22 degrees, that is, approximately 30 times 
better than COBE (7-10 degrees). Another major boost in data analysis is 
expected from the ESA satellite mission Planck, which is now scheduled to 
be launched on October 31, 2008; data releases for the public are expected 
in the following 3-5 years. Planck is planned to provide datasets of nearly 
5 x 10 10 observations, and this will allow to settle many open questions with 
CMB temperature data. New challenging questions are expected to arise at 
a faster and faster pace over the next decades; for instance, Planck will pro- 
vide high quality for so-called polarization data, which will set the agenda for 
the experiments to come. Polarization data can be viewed as tensor-valued, 
rather than scalar, observations — that is, what we observe are not measure- 
ments of a scalar quantity such as the temperature, but random quadratic 
forms. As such, this entails an entirely new field of statistical research, which 
is still in its infancy and will not be discussed in the present paper. 

Our aim here is to provide a review of statistical issues arising in CMB 
data analysis, with many examples of applications of statistical procedures 
to real data from the WMAP experiment. Some of the empirical results we 
provide are new, as detailed below. The plan of the paper is as follows: in 
Section 2 we review very briefly some background material on map-making, 
component separation and spectral representations for the CMB data sets. 
For brevity's sake, we do not provide many details other than the material 
which is essential for our following discussion. In Section 3 we are concerned 
with angular power spectrum estimation, and we discuss procedures to deal 
with relevant practical questions such as the presence of observational noise 
and/or missing observations. In Section 4 we present some tools to test 
for Gaussianity and/or isotropy of CMB radiation: we focus, in particular, 
on harmonic methods such as the bispectrum, techniques based on differ- 
ential geometry such as the local curvature, and spherical wavelets (with 
the so-called Spherical Mexican Hat approach). Concerning the latter, we 
stress that many other possible approaches to wavelets on the sphere exist, 
which have been successfully applied to various parts of cosmological and 
astrophysical research: nevertheless, the field is still extremely active and 
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very much open for research (in particular, the derivation of the stochastic 
properties of wavelets procedures is still at the very beginning) . Finally, we 
collect in the Appendix some background mathematical material which we 
considered necessary for a better understanding of our proposals. 

2. Some preliminary issues. 

2.1. Map-making and component separation. To understand more pre- 
cisely the nature of the statistical issues involved, we need to introduce 
some more formalization. As explained above, CMB can be viewed as the 
single realization of a random field on the surface of the sphere, that is, for 
each ieS 2 , T(x) is a random variable on a probability space. Observations 
are provided by means of electromagnetic detectors (so-called radiometers 
and/or bolometers) which measure fluxes of incoming radiations (i.e., pho- 
tons) on a range of different frequencies. For instance, the above mentioned 
WMAP experiment is based upon 16 detectors, centered at frequencies 40.7, 
60.8 and 93.5 GHz, which are labeled the Q, V and W band, respectively. 
The forthcoming ESA mission Planck will be based upon 70 channels rang- 
ing from 30 GHz to 857 GHz. As the satellites scan the sky, observations 
are collected as a vector time series, the number of observations being in 
the order of 10 9 for WMAP and 5 x 10 10 for Planck. A first issue then re- 
lates to the construction of spherical maps starting from the Time Ordered 
Data (TOD) provided by the satellite; this is the so-called map-making chal- 
lenge; see, for instance, Keihanen, Kurki-Suonio and Poutanen (2005) and 
De Gasperis et al. (2005). For brevity's sake, we shall provide only the ba- 
sic framework, and refer to the literature for more details. In short, we can 
assume that in each of the p channels we actually observe 

Oi(x) = T(x) + Fi(x) + Ni(x), i = l,...,p,xeS 2 ; 

here, T(-) denotes the CMB signal, Fi(x) denotes the so-called foreground 
emissions by galactic and extragalactic sources of noncosmological nature 
(for instance, galaxies, quasars, intergalactic dusts and others), and Ni(x) 
instrumental noise. The crucial point to be understood is that the depen- 
dence across the different frequency channels of CMB emission is known, and 
it is different from the pattern followed by other sources: this capital property 
makes component separation possible and allows the construction of filtered 
maps [see. e.g., Patanchon et al. (2005) and the references therein]. More 
precisely, a clear prediction from theoretical physics, confirmed to amazing 
accuracy from the very first experiments [Smoot et al. (1992)], is that CMB 
radiation should follow the Planckian curve of blackbody radiation, that is, 
radiation is distributed across frequencies z/j, i = 1, . . . ,p according to the 
function 

. . 8Trhv 3 1 

W R{v;x)= ^ e _ hv/kBT{x) _ r 
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Fig. 1. CMB radiation from W MAP data. 

where R{y\ x) denotes the emission at frequency v for the corresponding 
temperature T(x) (measured in Kelvin degrees), c is the speed of light in 
the vacuum (= 2.99798 x 10 8 m/s), h is Planck's constant (= 6.6261 x 1CT 27 
er g/s), and fee is Boltzmann's constant (= 1.3807 x 10~ 16 er g/K). In other 
words, the determination of T{x) is made possible by the inversion of (1): 
the blackbody pattern can be estimated due to the presence of multiple de- 
tectors and the fact that astrophysics! emissions of noncosmological nature 
are characterized by a different pattern of dependence across frequencies. In 
some regions, however, foreground emissions are so strong that component 
separation is still a difficult statistical problem; several groups of cosmolo- 
gists are active in this field and a unique consensus solution has not been 
delivered yet. Moreover, in some areas of the sky (e.g., the Galactic plane, 
i.e., the line of sight of the Milky Way) the problem is considered to be 
largely unsolvable, so that there are missing observations in CMB maps 
(these unobserved regions are becoming, however, smaller and smaller with 
more refined experiments). In Figure 1 we report a CMB map constructed 
from (the Q band of) WMAP data; the missing region around the galactic 
plane is immediately evident. 

Full-sky maps can be constructed by weighted linear interpolation across 
different channels, but they are not considered fully reliable for data analysis, 
especially at high frequencies; we report this so-called ILC (Internal Linear 
Combination) map in Figure 2, see Bennett et al. (2003) for more details on 
its construction. 
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Fig. 2. The so-called Internal Linear Combination map from WMAP data. 

There are several other statistically interesting issues involved with the 
reconstruction of the scalar value T(x) from the vector- valued observations 
{Oi(x), . . . ,O p (x)}; actually the real experimental set-up is more compli- 
cated (and interesting) than this, because each location is observed un- 
evenly, that is, the scanning strategy is such that some regions are more 
accurately measured than others. Also, the contaminating noise can have a 
time-dependent structure [there is indeed strong evidence for long memory 
behavior, see, e.g., Natoli et al. (2002)]; the possible existence of noise corre- 
lation across different channels will be discussed below. These experimental 
features have sparked in the cosmological literature a very lively statistical 
debate on filtering and image reconstruction. We shall come back to some 
of these points later. 

2.2. Isotropy and spectral representation. In the idealistic case of no ex- 
perimental noise and perfect map-making, we can focus on the random field 
{T(x)}, assuming that it is exactly observed at each location on the unit 
sphere S 2 . A crucial assumption on CMB radiation is its isotropic nature, 

that is, T(-) = To </(•), where = denotes equality in distribution (in the sense 
of random fields) and g € SO (3) is any element of the group of rotations in 
i? 3 . More explicitly, the joint law of CMB radiation is assumed invariant to 
any change of coordinate; the condition is viewed by the physicists as a real- 
ization of so-called Einstein's Cosmological Principle, that is, the statement 
that the Universe should "look the same" to an observer in any arbitrary 
location. In other words, we could impose isotropy by requiring that the 
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stochastic laws of CMB radiation are invariant with respect to the choice of 
coordinates. There is some (quite inconclusive) evidence from WMAP data 
that isotropy may fail, that is, some authors have suggested that data on 
CMB radiation may show some asymmetries which would be inconsistent 
with isotropy [see, e.g., Park (2004), Hansen et al. (2004)]. The existence of 
these asymmetries remains highly disputed, though, and it actually provides 
yet another intriguing area for statistical research. It is in fact hotly debated 
whether these asymmetries should be ascribed to experimental features or 
truly cosmological causes. From the theoretical point of view, cosmological 
models that would produce asymmetries do indeed exist, but they are highly 
nonstandard, ranging from global rotating solutions of Einstein's field equa- 
tions to unconventional topological structures for the whole Universe. Much 
more methodological and applied research is needed in this area, but the 
question will most probably remain unsolved at least until the first releases 
of Planck data are available in a few years' time. By now, it is fair to say that 
a vast majority of cosmologists is still sticking to the isotropy assumption, 
and this is what we shall do in the present paper. Some of the procedures we 
shall consider in Section 4 for testing non-Gaussianity, however, are known 
to have also power against nonisotropic behavior; see, for instance, the local 
curvature approach below. 

We shall hence focus on the statistical analysis of isotropic random fields. 
Throughout this paper we shall assume that the CMB random field is mean- 
square continuous, as it is always done in the CMB literature. Under the 
previous assumptions, the following spectral representation holds, in the 
mean square sense 

oo I 

(2) T{x)=Y J E a lmYl m (x) 

l=0m=-l 

(3) where ai m = T (x)Y i m {x) dx . 

Js 2 

Here, the bar denotes complex conjugation and {Yi m (-)} the spherical 
harmonics, which form an orthonormal system for L 2 functions on the 
sphere. Some explicit expressions for the spherical harmonics can be found 
in the Appendix: much more complete treatment can be found elsewhere; 
see Varshalovich, Moskalev and Khersonskii (1988). For I = m = 0, we have 
a oo = X52 T(x) dx, that is, the first coefficient is 4ir times the sample mean 
of the random field. This value can be subtracted from T(x), whence we can 
take the expansion to start from 1 = 1; indeed, in practice, in the cosmo- 
logical literature also the coefficients corresponding to / = 1 are discarded 
(the so-called dipole terms), as they have no cosmological meaning, but 
they simply reflect the absolute motion of the Earth with respect to the 
frame of reference with respect to which CMB radiation is at rest. For / > 2, 
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the triangular array {ai m (-)} represents zero-mean, complex-valued random 
coefficients, with variance E\ai m \ 2 = C\ > 0, the angular power spectrum of 



h. 



;m2 
mi ) 



the random field. The coefficients are uncorrected, Eai imi ai 2m2 = C^S^d', 
and, hence, in the Gaussian case they are independent [note, however, that 
o-im = (— l) m «i-m]- We have the identity 

{oo I \ 2 oo I 

Y^ Yl ai m Y lm (x) \ =J2Y1 E \ a lm\ 2 Yim(x) 
1=2 m=-l J 1=2 m=-l 

1 21+1 



J2Q J2 Y lm {x)=Y,Ci 



47T 
1=2 m=-l 1=2 

in view of a standard summation formula for spherical harmonics [Var- 
shalovich, Moskalev and Khersonskii (1988)]. It follows immediately that 
Ci(2l + 1) must be summable to ensure finite variance. The angular power 
spectrum in the Gaussian case provides a complete characterization of the 
dependence structure of the random field; to its estimation from CMB data 
we now turn our attention. 

3. Angular power spectrum estimation. 

3.1. Power spectrum estimation under idealistic circumstances. As noted 
before, having observed the random field T(x), the coefficients {ai m (-)} can 
be recovered by means of the inverse Fourier transform (3). In practice, 
with real data the integral is replaced by finite sums by means of (exact or 
approximate) cubature formulae, which are implemented in standard pack- 
ages for CMB data analysis such as HealPix or GLESP [see Gorski et al. 
(2005), Doroshkevich et al. (2005)]. The angular power spectrum can then 
be estimated by 

1 ' 

aim\ ■ 



(4) Q = 2FTT £ 



This simple estimator highlights a very important issue when dealing with 
CMB data. It is indeed readily seen that the estimator is consistent in the 
Gaussian case, as I — > oo; more precisely, 



ECi = Q 

1 

F, ■ 

Q 



(Q 1 2 1 



a l 



(2/ + 1) 

2TTT = °« 
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because af /Ci ~ Xi an d for m = 1, . . . , I, 2af m /Ci ~ i.i.d. x|> where Xn de- 
notes a standard chi-square random variable with n degrees of freedom. 
In the Gaussian case with fully observed maps, the issue of angular power 
spectrum estimation can thus be considered trivial, and indeed, the pre- 
vious expressions not only ensure consistency but they also provide exact 
confidence intervals: it is immediate to see that 

2 




2\a lrn \ z \ ~ Ci x xL+i- 

However, we must stress that these results rely heavily on the Gaussian 
assumption. Indeed, Baldi and Marinucci (2007) and Baldi, Marinucci and 
Varadarajan (2007) have shown that under isotropy the coefficients ai m can 
only be independent in the Gaussian case, despite the fact that they are 
always uncorrelated by construction: in other words, sampling independent, 
non-Gaussian random coefficients to generate maps according to (2) will 
always yield an anisotropic random field. The correlation structure of the 
coefficients {a/ m } is in general quite complicated, despite the fact that it can 
be very nicely characterized in terms of group representation properties for 
SO (3) [Marinucci and Peccati (2007)]. In view of this, to derive any asymp- 
totic result for G\ under non-Gaussianity is by no means trivial; indeed, even 
the possible consistency (as I — > oo) of the estimator (4) in non-Gaussian cir- 
cumstances is still an open issue for research. 

3.2. Dealing with instrumental noise. We shall now try to make our anal- 
ysis more realistic by considering the effect of noise and missing observa- 
tions. Starting from the former, we shall consider the case where we observe 
0(x) :=T(x) -\-N(x), N(x) denoting instrumental noise; for simplicity, we 
shall follow the cosmological literature, assuming N(x) to be also a zero 
mean, mean square continuous and isotropic random field on the sphere. 
Whereas the assumptions of zero-mean and mean square continuity are ba- 
sically immaterial, isotropy of the noise may need to be relaxed if the sky is 
unevenly observed. We shall also assume that T(x) and N(x) are indepen- 
dent. Performing the spherical harmonic transform, we obtain, in an obvious 
notation, 



aim= / {T(x)+N(x)}Y lm (x)dx=:af m + a^ n , 
Js 2 



which leads to 
1 



a 



21 + 1 



■ i i r i 

E l«L| 2 + E K| 2 + 2Re E "frrtfL 
.m=—l m=—l \m=—l 
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It is immediate to see that the resulting estimator is biased, ECi = Cf + Cj ; 
the variance is easily seen to be given by 

(5) VaX ^= 2{ %tf } - 

In the cosmological literature, the standard procedure to address this bias is 
to assume that the noise correlation structure can be derived by Monte Carlo 
simulations or instrumental calibration; under this assumption, it is possible 
to subtract the bias from C\ and obtain a correct estimator with variance 
(5). An obvious question is then to test whether the assumption that Cf 
is known does not introduce some spurious effect into the analysis (namely, 
some unaccounted bias). A proposal in this direction was put forward by 
Polenta et al. (2005). To understand this idea, we must get back to the 
multi-channel setting, where we observe 

O l {x):=T(x)+N l (x), i = l,...,p, 

which in the harmonic domain leads to 

T , Ni 

a i;lm ■— O-lra + a l m - 

Note that the temperature component of the random spherical harmonics 
coefficients does not depend on the observing channel. We assume that the 
noise is independent over channels, which is believed to be consistent with 
the actual experimental set-ups of current datasets. Testing noise correlation 
across different channels is yet another open challenge for research. For a 
given noise structure, an obvious estimator for C\ is 



Cf:=l±{Cu-CF}, 

(6) 

1 
:= o; I i 2-^i l a M m l • 

m=— I 

The estimator Cf- is known in the literature as the auto-power spectrum. 
Simple computations yield [Polenta et al. (2005)] 

EC l =Ci, 

Of course, the natural question that arises at this stage is the possible ex- 
istence of misspecification, that is, some errors in the bias-correction term 
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C h *. A solution for this issue was proposed by Polenta et al. (2005). The 
idea is to focus on the cross-power spectrum estimator 



'■iCP 



2 £i A / 



r " ' !=1 j=!+l \ m=—l / 

The underlying rationale for Cp is easy to gather: under the assumption 
that noise is independent across a different channel, the estimator is unbi- 
ased, regardless of the value of the C l i . More precisely, 



2 P " X p 



TyEE^ 



p(p - 1) • 



Ci. 



Similar manipulations yield 

o ( or< p -, p-i p 

2/ + 11 ' P 2 ^' P*(P-l) 2 £f,±ii ' ' 



Merely for notational simplicity, we also assume that the noise variance is 
constant across detectors. It is then readily seen that 

More explicitly, the auto-power spectrum estimator is more efficient that the 
cross-power spectrum; however, the latter is robust to noise misspecification. 
This is the classical setting which makes the implementation of a Hausman- 
type test for misspecification feasible [Hausman (1978)]. Indeed, it is possible 
to consider the statistic 

ff, = [Var{Cf * -CftrV'lCf -C?}, 

Under the null of exact bias correction, it is readily seen that H[ — >^ N(0, 1), 
as I — > oo. On the other hand, in the presence of misspecification, that is, 
when the actual noise variance is equal to C ; l + 5 for some i, 5 > 0, then 
we expect EHi to diverge with rate y/l5 as I — > oo. 

It is also possible to consider a functional form of the same test, focusing 
on 

1 [Lr] 

B L {r):=^=Y. H ^ r€[0,l}. 
vL l=1 
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It is standard to show that Bi(r) converges weakly to a standard Brown- 
ian motion, as L — ► oo. A test for noise misspecification can then be con- 
structed along the lines of standard Kolmogorov-Smirnov or Cramer- Von 
Mises statistics. We refer again to Polenta et al. (2005) for a much more 
detailed discussion and an extensive simulation study. 

The methods discussed above rely on a basic identification assumption, 
that is, the condition that instrumental noise be independent across dif- 
ferent channels. This is an assumption which is commonly entertained in 
the cosmological literature; suitable statistical issues to test its validity are 
still lacking and represent an open issue for research. A more challenging 
research task was mentioned before: the previous discussion was entirely 
led under the assumption that the CMB field (and thus the corresponding 
spherical harmonics coefficients) are Gaussian. It is very important to stress 
that relaxing this assumption has much deeper consequences here than it 
is usually the case in statistical inference. Indeed, it follows from results 
in Baldi and Marinucci (2007) that if the field is isotropic, the coefficients 
(ai m ) cannot be independent unless they are Gaussian. It follows that even 
the simple consistency (as / — > oo) of the estimator C\ remains an open issue 
to address, in general non-Gaussian circumstances. We shall not go further 
into this issue here, but we rather focus on another important feature of 
realistic datasets: the presence of unobserved regions, which make the exact 
evaluation of the inverse Fourier transform (3) unfeasible. 

3.3. Missing observations. The presence of missing observations, that is, 
regions of the sky where the CMB is deeply contaminated by astrophysical 
foregrounds, posits serious challenges to angular power spectrum estimation. 
The first consequence is that the sample spherical harmonics coefficients 

M 



a im= / T(x)Y lm (x)dx, 

JS 2 /M 

lose their uncorrelation properties (here, M denotes the unobserved region 
and, for notational simplicity, we came back to the case of a single detector 
with no instrumental noise). Indeed, we have 



Eaf? mi a% m2 = E{ ( / T(x)Y hmi (x) dx^j (^ ^ T(y)Y hrri2 (y) dy^j j 
V V Eai m ai' m > ( / Yi m (x)Y hmi (x) o 

rz rz \Js 2 /M 



(7) 



hm 1 l 2 m 2 JS 2 /M 

x ( / Y Vm i(y)Yi 2m2 (y)dy 

\JS 2 /M 

2J CiWi m i imi Wi m i 2m2 , 

lin 
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where 



Wi m i imi := Yi m (x)Y hmi (x)dx. 

JS 2 /M 

In case the spherical random field is fully observed, then M = (the empty 
set) and by standard orthonormality properties of the spherical harmonics 
Yi m , we obtain Wi m i imi = ^d and > therefore, Ea hmi ai 2rri2 = Ci6\*5™£. I n 
the presence of missing observations, the random coefficients are no longer 
uncorrelated neither over I nor over m. In the physical literature the val- 
ues of {Wi m i imi }i imi i 2m2 are computed numerically, exploiting the a priori 
knowledge on the geometry of the unobserved regions; the resulting coupling 
matrices can then be used to deconvolve the estimated values Ci, a proce- 
dure which has become extremely popular under the name of MASTER [see 
Hivon et al. (2002) for details]. In practice, it is not possible to identify by 
this method the value of the angular power spectrum at every single multi- 
pole I; it is then customary to proceed with binning techniques, where the 
values of C\ at nearby frequencies are averaged and only these smoothed val- 
ues are actually estimated. Plots for the estimates of the C\ derived along 
these lines can be found, for instance, on the web site of WMAP; a compar- 
ison with angular power spectrum estimate from several other experiments 
(based upon smaller patches of the observed sky) is also entertained. 

The previous procedures can be computationally extremely demanding 
and we would like here to introduce an alternative strategy, which was 
basically put forward in Baldi et al. (2006). The idea is to implement 
power spectrum estimation by means of new kinds of spherical 
wavelets, called needlets [see also Narcowich, Petrushev and Ward (2006a), 
Narcowich, Petrushev and Ward (2006b), Marinucci et al. (2008) and 
Baldi et al. (2007)]. Needlets can be described as a convolution of the spher- 
ical harmonics basis by means of a suitable kernel function &(•); more pre- 
cisely, the general element of the needlet frame can be written down as 

I BJ+1 ( I \ l - 

J=BJ-i ^ ' m=-l 

where {£■,&} denotes a set of grid points on the sphere, B > 1 is a bandwidth 
parameter, &(•) is compactly supported and an infinitely differentiable func- 
tion which satisfies the partition- of-unity property, that is, 

(8) E &2 (^j) El foralU>l, 

3 

and {\jk,£,jk} (the cubature points and cubature weights) can be chosen in 
such a way that 

J2 Y hmi((,jk)Yi 2m2 ((,jk)>>jk= / Y hmi (x)Yi 2m2 (x)dx = 5\l5™l. 
k Js2 
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More details on this construction and its underlying rationale can be found 
in Baldi et al. (2006) and are not reported here for brevity's sake; see also 
Kerkyacharian et al. (2007) and Guilloux, Fay and Cardoso (2007) for fur- 
ther work in this area. The corresponding random needlets coefficients are 
provided by the analysis formula 

r I — BJ+1 l ( I \ - 

Pjk= / T(x)ip jk (x)dx = JXjk V V b — )ai m Y lm (£j k ), 

J V l=B^m=-l ^ WJ 

whereas the synthesis expression is given as 

BJ+1 h ( h \ ( h \ 

J2^^jk(x)=Y^ E E b (-E]) b (w) ahmiYlimi ^ 

j,k j i 1= Bi-i mi=-Zi V-° / \-° / 

BJ + 1 h 

X Yl Yl Y Y hrni(£jk)Y 'l 2 m 2 {£jk)Xjk 

l 2 =Bi- 1 m2=-h k 

BJ+ 1 I 



j l 1 =Bi- 1 m 1 =~l 1 V / V / 



J h=B3 

BJ + 1 h 

x E E %%£ 

l 2 =Bi- 1 rn2=-h 
oo / 

= E E a imYi m {x)=T(x), 

1=1 m=—l 

using (8). For our purposes, it is sufficient to recall the main properties of 
the needlets construction: 

• needlets enjoy excellent localization properties in the real domain, each 
ipjk(x) being quasi-exponentially localized around its center £j&. As such, 
needlets coefficients have been shown to be minimally influenced by the 
presence of missing observations. 

• the needlets system is compactly supported in the harmonic domain; as 
such, the random needlets coefficients are uncorrelated for j — j' > 2. Much 
more surprisingly, the random needlets coefficients are asymptotically un- 
correlated for any fixed angular distance, as the frequency j diverges to 
infinity. This property implies that (in the Gaussian case) it is possible to 
derive a growing array of asymptotically i.i.d. observations out of a single 
realization of an isotropic random field. This opens the way to a plethora 
of statistical procedures. 
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In particular, it is possible to suggest the estimator 

r BJ+1 l i — / 1 \ i 2 

r ^=E^=E E E v^ 6 (hj )«'m*w&k) 

= E E b [w) b [w) aiimiSl2ma \ A J fe E y ^(^' fe ) y ^fe fe ) r 
= E a^cica + i), 

for which it is simple to show that 

BJ+1 / I \ 

(9) sr i= e h (-hm 2l+l )- 

Equation (9) shows that Tj provides an unbiased estimator for a smoothed 
version of the angular power spectrum; the advantage with respect to the 
standard procedure is that not only unbiasedness, but even uncorrelation 
over different scales is asymptotically conserved in the presence of missing 
observations, making the implementation of confidence intervals and test- 
ing procedures viable [see again Baldi et al. (2006) for details]. Also, even 
in the presence of a masked region, the summands {/3jk} are still asymp- 
totically independent (over k) as j — > oo, whereas we have seen in (7) that 
this is not the case for the random coefficients {aim}. The price for such 
robustness properties is clearly connected to the smoothing, that is, in the 
presence of missing observations it turns out to be unfeasible to estimate 
each angular power spectrum mode C\ by itself, and one must stick to a 
slightly less ambitious goal, that is, the estimation of joint values averaged 
over some subset of frequencies (chosen by the data analyst). There is, of 
course, a standard trade-off in the choice of the bandwidth parameter B: 
values closer to unity entail a much better resolution, but this brings about 
worse localization properties on the sphere and therefore a possibly higher 
contamination from spurious observations; on the other hand, higher values 
of B yield more robust, but less informative estimates. 

Spherical wavelets in general, and needlets in particular, allow for many 
statistical applications, which go much beyond angular power spectrum esti- 
mation. One example is the analysis of cross-correlation between CMB and 
Large Scale Structure (LSS) maps; this is a key prediction of many cosmolog- 
ical models entailing some form of dark energy and has been implemented on 
real data by Pietrobon, Balbi and Marinucci (2006). Other applications may 
include testing for non-Gaussianity and isotropy, bootstrap/subsampling 
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evaluation of confidence intervals for CMB statistics [Baldi et al. (2007)], 
component separation and many others. Given such a wide array of appli- 
cations, we stress the need for a more careful analysis of their theoretical 
underpinnings, with special reference to the effect of the Gaussianity as- 
sumption on our conclusions. This and many other related issues are left as 
topics for further research. 

3.4. Parameter estimation. In this paper we shall neglect almost com- 
pletely another crucial issue in CMB data analysis, which is very tightly 
coupled to the estimation of the angular power spectrum, that is, cosmo- 
logical parameter estimation. More precisely, the theoretical angular power 
spectrum can be written as a function of a number of cosmological param- 
eters, such as the baryon, matter and dark energy densities Qfe,f2 m ,QA, the 
optical depth r, the spectral index n s , the Hubble constant Hq and others; 
of course, the numbers of parameters to be estimated varies across differ- 
ent cosmological models, typically ranging from 6 to 16; see again Dodelson 
(2003) for more details. There are no known closed-form expressions yielding 
the theoretical angular power spectrum C\ as an explicit function of these 
parameters (which we write for brevity as #); however, there are indeed 
very fast numerical routines which solve the associated partial differential 
equations and provide as an output C\ after a specific value of $ has been 
supplied [see Seljak and Zaldarriaga (1996)]. 

Once the set of estimated values C\ has also been derived, there are ba- 
sically two approaches that have been implemented to obtain estimates for 
the set of parameters, namely, some form of minimum distance estimators, 
where the parameters are calibrated to minimize a weighted distance be- 
tween C;(i?) and Ci, and approximate maximum likelihood methods, where 
suitable approximations for the likelihood functions are derived and the esti- 
mates are consequently derived. In practice, both methods are implemented 
by means of a heavy use of numerical techniques (especially MCMC), and 
a lively debate is growing on the construction of the most efficient algo- 
rithms. Likewise, an extensive discussion is growing on the construction of 
confidence intervals for the parameters, where fundamental issues such as 
the differences between Bayesian and frequentist viewpoints are often called 
upon (the distinction between these two approaches is not perceived in the 
cosmological community in the same manner as in the statistical one; just 
to give an example, maximum likelihood estimates are nearly unanimously 
labeled a Bayesian procedure in the CMB literature). 

For brevity's sake, we are unable to go deeper into these issues, which 
are still quite far from a satisfactory solution. We refer, for instance, to 
Hamann and Wong (2008) and the references therein for more discussion 
and recent proposals in this area. 
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4. Testing for non-Gaussianity. Among the several statistical issues which 
arise in connection to the analysis of Cosmic Microwave Background radia- 
tion, a lot of attention has been drawn by non-Gaussianity tests. These tests 
have several motivations. The first is connected to the need for a statistical 
validation for the predictions of the so-called inflationary scenario, which 
is currently the leading incumbent as a standard model for the Big Bang 
dynamics; see Dodelson (2003) for discussions and explanations. Under this 
labeling, there exist an enormous variety of different physical models, which 
in a vast majority of circumstances lead to expressions such as 

(10) T(x) = T G (x) + f N L{Tl(x) - ET^(x)}, 

where T(x) denotes as before CMB, Tq{x) is an underlying Gaussian field, 
/jvl is a nonlinearity parameter and the unit of measurements are such that 
the non-Gaussian part Tq(x) — ETq(x) is 10~ 4 /10 -5 times smaller than 
Tg(x). (10) should be viewed as a strong simplification, for several reasons: 
in particular, we are considering exclusively the primordial dynamics, thus 
neglecting later interactions through the gravitational potential; also, we are 
ruling out more complicated models, where higher order terms or multiple 
subordinating fields may be present; and, of course, we are neglecting a whole 
plethora of observational issues, where possible non-Gaussianities may be 
formed by secondary effects, such as the interactions of incoming photons 
at more recent epochs. Despite all these simplifying conditions, (10) does 
provide an extremely good guidance for features to be expected and, indeed, 
it makes up a benchmark model against which many procedures have been 
tested in the last few years. In particular, considerable attention has been 
drawn by the possibility to constrain the value of /jvxi as this depends on 
constants from fundamental physics [Bartolo et al. (2004)] and as such it 
allows to probe many features of cosmological models. 

Among several statistical procedures which have been proposed in the 
literature, we shall focus on three main families, namely, tests based upon 
the bispectrum, tests based upon geometric features of Gaussian random 
fields (local curvature) and tests based upon spherical wavelets (in this case, 
so-called Spherical Mexican Hat Wavelets). 

4.1. The angular bispectrum. It is obvious that, under Gaussianity, the 
sequence {ai m } , m = 0, . . . , I makes up an array of independent Gaussian ran- 
dom variables (complex- valued for m ^ 0), so that a natural first option for a 
test of Gaussianity is to consider their sample skewness ai imi ai 2m2 ai 3m . i and 
check whether it is significantly different from zero. This simple idea is made 
much more sophisticated by the necessity to impose rotational invariance on 
the sample coefficients. Such invariance can be imposed by demanding that 
the probability law of the CMB field be invariant with respect to the action 
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of the rotation group. More formally, let g £ SO (3) be any element of the 
rotation group in M 3 ; the assumption of isotropy can then be written as 

T(x) = T(gx) for all x G S , 



whereas in terms of the spectral representation, we have 

oo I OO I 

(11) E E ai m Y [m (x)=Y E a lm Y lm(gx). 

1=0 m=—l 1=0 m=—l 

As explained, for instance, in Hu (2001) and Marinucci and Peccati (2007), 
from (11) it follows that the bispectrum of a rotational invariant random 
field must take the form 

(h l 2 h\(h h h\ h 
ua hmi a l2m2 ai 3rn3 - y Q Q Q j ^ ^ m3 J°hhh, 

where 6j 1 ; 2 ; 3 (the reduced bispectrum) conveys the physical information and 
does not depend on 7711,7712,7713. The Wigner's 3j symbols appearing on 
the right-hand side are discussed in the Appendix; many more details can 
be found, for instance, in Varshalovich, Moskalev and Khersonskii (1988), 
Marinucci (2006) and Marinucci (2008), whereas generalization to higher 
order cumulant spectra are described in Marinucci and Peccati (2008). A 
feasible, rotationally invariant estimator for the (normalized) bispectrum is 
provided by 

Q>limi 0; 2 m2 ^'h rn 3 



4^3 = (-i) (/l+/2+/3)/2 E (t 

mim2m,3 ^ 


h 

77l 2 


h 
m 3 


and its studentized version is of course 






4^ = (-i) (h+ ' 2+/a)/2 E (1 


h 

m 2 


h 
m 3 



^C h Ci 2 Ci 3 



"limi O j l2ni2 ^£37713 



The sample bispectrum is discussed, for instance, by Hu (2001); asymptotic 
properties are provided by Marinucci (2006) and Marinucci (2008), where the 
phase factor (— l)('i+' 2 +'3)/ 2 j s a i so introduced. In particular, it can be shown 

that the sequences {iw 2 Z 3 } an d {Ihhh} converge to Gaussian independent 
random variables in the high frequency limit where min(Ji,l2,h) T °°- The 
limiting behavior of the bispectrum ordinates, however, is perhaps not the 
most significant instrument for the implementation of statistical procedures. 
More precisely, it seems more promising to combine the different ordinates 
into a single statistic, by means of the integrated bispectra 

[Lr\ ( K ^ [Lr] ^ 

JlL(r)= El -7^Y,^l+k,hh h J2L(r) =^2Tiu. 

h=1 {VK k=1 J l=1 
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Convergence to Brownian motion for both these statistics was established 
in Marinucci (2006). The underlying rationale can be briefly motivated as 
follows: in both cases we combine several different ordinates into a single 
functional statistic, capable of keeping track of the frequency location for 
possible deviations from Gaussianity. The different combination of multi- 
poles in Jil(-), J2l(-) corresponds to the two well-known classes of squeezed 
and equilateral configurations, as discussed again by Marinucci (2006), 
Babich, Creminelli and Zaldarriaga (2004) and many others. It is also pos- 
sible to provide some results on the asymptotic behavior of these statistics 
under non-Gaussian circumstances; in particular, results in Marinucci (2006) 
suggest that J\l will provide consistent testing procedures (as L — > oo) under 
model (10), whereas tests based upon J2L will have asymptotically negligible 
power, for all values of /jvl- These theoretical findings have been validated 
by Monte Carlo simulations in Cabella et al. (2006); the integrated bispec- 
trum has also been shown to compare favorably with alternative statistical 
procedures in some internal statistical challenges within the Planck collab- 
oration. 

In Figure 3 we report the results obtained by implementing Jix(r) on 
the data from the (2003) and (2007) WMAP data releases. We stress that 
the simulations are calibrated in a realistic experimental setting, that is, 
they do take into account features such as the presence of noise and missing 
observations. More precisely, we used 1000 simulated maps of CMB signal 
plus noise; we took into account the modulation of noise on the maps given 
by WMAP scanning strategy, the presence of a masked region to avoid the 
emission from the Milky Way and point sources, and we considered the 
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Fig. 3. The behavior of Jii(r) on WMAP data. 
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optical transfer function of the telescopes. To comply with the cosmological 
literature, the shaded region represents the 68% confidence interval (la) 
as evaluated by means of Monte Carlo simulations for various values of 
r € [0, 1]: we fixed L = 500 because WMAP data allow a reliable coverage 
up to this multipoles; see Bennett et al. (2003) and Hinshaw et al. (2008) 
for more technical details on the experiment (it should be noted that in the 
x-axis we report rL). The dotted and dashed lines represent Monte Carlo 
expected values for our statistics with f^L = ±100, . . . , 500, respectively. It 
is possible to check that the boundary value of Jnl to ensure detection is in 
the order of 200 or larger, that is, with a signal to noise ratio in the order of 
a few percentage points. This is indeed confirmed by a more detailed study 
in Cabella et al. (2006). Finally, triangles (2003 dataset) and squares (2007 
dataset) represent the evaluation of the statistic on real data, on the basis 
of the previously mentioned WMAP releases. It is clear that the evidence 
for non-Gaussianity is rather weak, and, indeed, the statistics get closer to 
zero as the observations increase. We must stress, however, that the level 
of non-Gaussianity favored by theorists is well below 100, and this is still 
consistent with observations at the current resolution. Note that the signal 
to noise ratio for the non-Gaussian signal is in the order of /wx/lO 4 , so that 
these values are really difficult to detect. 

Very recently, in Yadav and Wandelt (2007) it has been claimed a detec- 
tion of a nonzero /nl (— 80) by means of a modified bispectrum estimator, 
which is constructed to take into account the presence of noise and miss- 
ing observations, at the same time keeping computational costs at a feasible 
level. This proposal is indeed very interesting; the results, however, are quite 
close to the boundary level and as such they must probably be considered not 
conclusive. The general consensus in the community seems to be that new 
releases of data from more sophisticated experiments such as Planck, and 
possibly more efficient statistical procedures yet to be devised, will indeed be 
necessary to settle the question on the possible existence of non-Gaussianity 
in CMB. It should be stressed, in particular, that the bispectrum requires 
the evaluation of the inverse Fourier transform (3), and as such it is known 
to be severely affected by the presence of missing observations (there is some 
evidence that the detection level could reach _f/vx — 10 or lower for fully ob- 
served sky maps). Improving the performance of the bispectrum for partial 
sky coverage is a priority of current research in view of the forthcoming 
satellite data: for instance, in Lan and Marinucci (2008) the bispectrum ap- 
proach is combined with the needlets construction described in the previous 
section. Rather than considering these further developments in the bispec- 
trum literature, we move to other methods which have a local nature, and 
are thus expected to be more robust in the presence of missing data. 
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4.2. Local curvature. The next approach to constraining possible non- 
Gaussianities is based upon an analysis of the local features of Gaussian 
random fields. This is an issue which has a long tradition in probability, as 
summarized, for instance, in the recent book by Adler and Taylor (2007). 
In this respect, many efforts have focused on excursion sets and other pro- 
cedures from convex geometry (the so-called Minkowski functionals) ; these 
ideas have found many fruitful applications in CMB data analysis; see, for 
instance, Hikage et al. (2008) and the references therein. 

However, our approach here will be different, and more directly rooted to 
differential geometry; we collect in the Appendix some background material 
to understand better the notation and the approach. The intuition can be 
explained as follows. At any point x £ S 2 on the sphere, it is possible to 
investigate the local curvature of the random field T by focusing on its 
covariant Hessian matrix; in particular, we can study whether this Hessian 
defines a positive definite bilinear form (in which case we will label x as 
lake point), a negative definite form (in which case we will label x as hill 
point) or neither of the two (in which case we will label x as saddle point). 
This approach was proposed by Dore', Colombi and Bouchet (2003), in the 
standard Euclidean circumstances, and then applied to the spherical case 
by Hansen et al. (2004) and Cabella et al. (2005). Here, a crucial step is to 
ensure that the Hessian is defined in such a way to have an intrinsic meaning, 
that is, the geometric characterization of the points must be independent 
from the choice of coordinates. The appropriate instrument for this point is 
the notion of covariant derivative, which we recall briefly in the Appendix. 
We are finally in the position to evaluate the covariant Hessian of any random 
field on the sphere, which is provided by 

\T.a v / sin$ T :w> /sin 2 ?? 
(12) ] ' V 

T## (T& v - cot i?T j¥ ,)/ sin •& 

{Ttfy — coti?T i¥ ,)/sini? (T !tpip + sin?9cos??T,9)/sin 2 $ 



where for a,b = , &,<p we have 
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Explicit expressions for partial derivatives of the spherical harmonics can 
be found in Varshalovich, Moskalev and Khersonskii (1988); we have, for 
instance, 

d 

— Yim^, if) = imYirnitf, if), 

dip 



^Y lm {<&,ip) = l -^l{l + 1) - m(m + l)Yi, m+ i(tf, ?>'* 



22 P. CABELLA AND D. MARINUCCI 

- l^l(l + i)- m(m - l)Y|, ro _i(0, <p)e*. 

In particular, the eigenvalues of the Hessian matrix H are intrinsic quanti- 
ties and do not depend on the choice of coordinates; hence, points can be 
classified as hills, lakes or saddles by checking whether both eigenvalues are 
positive, both negative or with opposite sign, respectively. The next step in 
the procedure is to focus on thresholded random fields, that is, to consider 
only those values where T > v, for some level v which is taken between ±3cr 
from zero (a denoting as usually the standard deviation of T). The relative 
proportion of any two of the curvature typologies can then be evaluated as 
a function of a threshold value v\ that is, if we consider hills h(v) and lakes 
l(v) we get 

#{xi : X 1 (H(x i )),X 2 (H(x i )) > 0, T(xi) > v] 



h(v) :-- 



l(v) :- 



#{xi-.T(xi)>v} 
#{xi : AiOff(xi)), X 2 (H(xi)) < 0,T( Xi ) > u] 



#{xi:T(xi) >v} 

where \i(H(xi)),\2(H(xi)) denote the two eigenvalues of H(-) at the lo- 
cation Xi, {xi} denoting any discretization of the sphere as provided, for 
instance, by HealPix. The same procedures can then be evaluated on a grid 
of different threshold values Uj, j = l,...,q, and this leads to normalized 
statistics 

i' (l .)- l M h'(v)- h{Vj) 

It must be stressed that Dore', Colombi and Bouchet (2003) provided some 
analytic results for El(vj),Eh(vj) in the standard Euclidean case; as these 
procedures depend only on local features, these analytic results provide ex- 
cellent approximations even in the spherical case, as shown in Cabella et al. 
(2005). On the other hand, currently there is no rigorous result on the asymp- 
totic distribution of such statistics, which must hence be calibrated by sim- 
ulations. 

In Figures 4 and 5 we report the la confidence regions for the hills and 
lakes densities at various thresholds, evaluated as before by simulations on 
1000 Gaussian random fields which mimics the expected behavior of CMB 
radiations [see again Cabella et al. (2005) for details]. As in the previous 
subsection, the dotted and dashed lines represent Monte Carlo expected 
values for the values of our statistics for f^i = ±100, . . . , 500, respectively. 
We also report our estimates based on the WMAP 2003 (crosses) and 2007 
(squares) data releases. The general conclusions seem rather close to what 
we found for the bispectrum: the evidence for non-Gaussianity is apparently 
weak. On the other hand, the non-Gaussian simulations seem to suggest 
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Fig. 4. "Hills" density on WMAP data. 

that the power here may be slightly weaker than for the bispectrum, and in 
any case insufficient to detect values of /jvl smaller than 100, as predicted 
by the theorists. Again, the new data releases from Planck are mandatory 
to reach firmer conclusions in this area. 

As a final remark, we stress that local curvature methods entail a possi- 
bility which is ruled out by the bispectrum: as the methods are local, they 
can be used to test for isotropy, for instance, comparing the behavior of 
the local curvature on different hemispheres of the CMB sky. This is in- 
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deed the approach pursued by Hansen et al. (2004), where the results on 
WMAP data are compared with Monte Carlo simulations, presenting some 
boundary evidence that the assumption of isotropy may fail. The literature 
on the possible existence of these asymmetries has grown enormously in the 
last 4 years, but no consensus has been reached. As mentioned earlier, the 
existence of asymmetries would have very profound consequences on cosmo- 
logical theory, and again, much progress in this area is awaited in the next 
decade. 

4.3. Spherical Mexican Hat Wavelets. Many different proposals have been 
put forward for the implementation of spherical wavelets systems on the 
sphere: for the approaches which are more directly connected to the cos- 
mological community, see, for instance, Antoine and Vandergheynst (2007), 
McEwen et al. (2007), Wiaux, McEwen and Vielva (2007) and the refer- 
ences therein. The construction by Antoine and Vandergheynst (1999), later 
developed by Wiaux, Jacques and Vandergheynst (2005) has become espe- 
cially popular; application to CMB data with the aim of testing for possible 
non-Gaussianity is due to Cruz et al. (2007) and Cruz et al. (2006). We shall 
focus here on a version of the same approach, that is, the so-called Spherical 
Mexican Hat Wavelets (hereafter SMHW) . The idea of the construction can 
be explained as follows: in general, a wavelet system on M can be character- 
ized by means of dilations and translations of a fundamental function (the 
mother wavelet). On the sphere, the idea is to replace the translations by 
rotations, that is, elements of the group SO (3). To implement dilations, we 
note that around the North Pole the latter can be implemented by consider- 
ing usual dilations in the tangent plane, which are lifted on the sphere by in- 
verse stereographic projections from the South Pole. More precisely, after the 
identification of the tangent plane at the North Pole with the complex line 
C, the projection of a point u> = ($,<£>) is provided by $(o;) =: C = 2e l{f> taxi -^ ■ 
So a stereographic dilation D a :S 2 — ► S 2 reads Z) a (#,</?) = (fia^), where 
d a :tan% = a tan | , that is, # a := arctan{2atan -^}. 

More explicitly, the procedure to implement SMHW can be described as 
follows. In R 2 , the continuous Mexican Hat Wavelet can be written as 



V(x,R): 



exp(-|x| 2 /2 J R), 



<2itR 
which satisfies the standard compensation and admissibility conditions 

Ir 2 Jr 2 

the hat denoting Fourier transform. For a given scale R and location x £ S 2 , 
the definition of the (continuous) Spherical Mexican Hat Wavelet transform 
can then be entertained in three steps: 



/ ^(x,R)dx = 0, f iJfe^l dx=:Cy<oo, 

Jr 2 Jr 2 x 
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• A change of coordinates is performed, to rotate x into the North Pole. 

• A stereographic projection on the tangent plane is implemented, so that 
y := 2tan( : |), d denoting as usual the radial distance from the North Pole. 
We then implement the MHW on y. 

• A rotation is performed to transform back the wavelet coefficients to the 
original location. 

It should be noted how the same formalism can be fully justified with- 
out the need for stereographic projections, and resorting instead to group 
representation theory; we refer to Antoine and Vandergheynst (1999) and 
Wiaux, Jacques and Vandergheynst (2005) for more details. Following this 
route, the SMHW basis that we implement is given at the North Pole by 



*(s,i2):= 



2ttN(R) 



1 + 



.r 



2 n 2 



X \ 

r) 



2-1 



exp(-x 2 /2R), 



with corresponding random coefficients 



w(R):= / T(x)V(x,R)dx, 
Js 2 

where N(R) is a normalizing constant and x denotes the polar angle obtained 
with the stereographic projection; see also Vielva et al. (2003), Cruz et al. 
(2007) and Cruz et al. (2006). SMHW do not represent a tight frame, so 
no exact reconstruction formula is available. Their stochastic properties are 
currently under investigation to establish whether their random coefficients 
are asymptotically uncorrelated, as it was the case for needlets. On the other 
hand, SMHW do enjoy very good localization properties [see Marinucci et al. 
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Fig. 6. The behavior of SMHW on WMAP data. 
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(2008) for a comparison with needlets], they are simple to implement and 
they have been very widely used for Gaussianity and isotropy testing. For 
completeness, we thus implement the resulting coefficients as a test of non- 
Gaussianity, to be compared with our previous findings. In particular, we 
considered skewness and kurtosis for the SMHW random coefficients, which 
we calibrated by means of Gaussian simulations. For the skewness, the result 
are reported in Figure 6; for kurtosis, they are not reported, as the resulting 
power properties where worse. The shaded area, the dotted and dashed 
lines, the crosses (2003) and the squares (2007) have the same meaning as 
before: it should be noted that in the x-axis is reported the scale factor 
R, so that when moving from left to right we are approaching larger scales 
(i.e., smaller frequencies, the opposite than for the bispectrum, which is a 
frequency-domain statistic) . 

As before, the evidence for non-Gaussianity is weak; worse than that, 
here simulations suggest that much larger values of Jnl would be needed to 
ensure detection. It seems thus that this class of methods cannot outper- 
form procedures such as the bispectrum when looking for non-Gaussianity. 
It must be recalled, however, that wavelets do enjoy a very important ad- 
vantage on pure harmonic methods: indeed, their localization properties in 
the real domain allow the detection of unexpected features which may signal 
anisotropic behavior. A striking example of this possibility is provided by 
Cruz et al. (2007) and Cruz et al. (2006), where a form of SMHW has been 
used to detect a Cold Spot in CMB radiation. The existence and possible 
explanations for such features are again very widely debated — there seems 
to be a tight relationship with the evidence on asymmetries which was men- 
tioned earlier [Park (2004) and Hansen et al. (2004)]. This is one more area 
where new statistical challenges will take place on Planck data, and the most 
suitable forms of spherical wavelets will certainly provide valuable contribu- 
tions. 

APPENDIX 

A.l. Isotropy and Wigner's coefficients. A proper derivation of the spher- 
ical bispectrum expression would require a considerable effort and some non- 
trivial background on group representation theory. We report here just the 
basic facts, and refer to the literature for a more detailed discussion [see, 
e.g., Varshalovich, Moskalev and Khersonskii (1988), Vilenkin and Klymik 
(1991), Hu (2001) and Marinucci (2006, 2008)]. 

The spherical harmonics are defined by 



111 + 1 (I — m)\ 

Yi m (9,ip) = «/ — — -P; m (cos6>)exp(imy?) for m > 0, 

y 47r (l + m)\ 

Y lm (9,<p) = (-l) m Y l]m] (9,<p) for m < 0, 
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where Pi m (cos9) denotes the associated Legendre polynomial of degree l,m, 
that is, 

p lm { x ) = (-i)-(i - x ^— Pl ( X ), Pl { x ) = w —^ - iy, 

m = 0,l,2,...,l, Z = 1,2,3, 

Here, < $ < ir and < ■$ < 2tt denote the usual spherical coordinates on 
the sphere; more explicit expressions for YJ m (-) are given below [see (14)]. 
The spherical harmonics provide an orthonormal system for the space of 
square integrable functions L 2 (S 2 ) with the uniform measure. Now let g £ 
5*0(3) be an arbitrary rotation of M 3 (i.e., a change of coordinate). It is a 
well-known fact in spherical geometry that the action of the rotation group 
can be parametrized by the three Euler angles g = (a,(3,j), < a < 2tt, 
0</3<7r, 0<7< 2-7T. The action of SO(3) on the spherical harmonics can 
instead be expressed as 

i 
(13) Yi m (gx)= ^2 Dl m'm(9)Yi m >(x), 

m'=— I 

where {D l m / m (-)} are the so-called Wigner's D matrices, which provide an 
irreducible representation of the group of rotations SO (3). In coordinates, 
an explicit expression for the elements of the D-matrices is provided by 

D l (a, (3, 7 ) = {D l m , m (a, (3, 7 ) V, m =-z,...,z = {e- %m ' a d l m , m ^)e m ^ m ^ m= . K .. h 
where 

d l m M = (-l) l - n [(l + m)\(l - m )\(l + n)\(l- n)!] 1 / 2 

^ k (cos (^/2)) rn+ra+2fc (sin(^/2)) 2; - m - n - 2fc 
X ^ ' k\(l -m- k)\(l -n- k)\{m + n + A;)! ' 

and the sum runs over all k such that the factorials are nonnegative; see 
Varshalovich, Moskalev and Khersonskii (1988) and Vilenkin and Klymik (1991) 
for a huge collection of alternative expressions. The proof of (13) is based 
upon group representation theory techniques and we do not provide it here; 
we simply recall that the elements of D l (a,(3,^f) are related to the spherical 
harmonics by the relationship 



(14) D l 0m (a, A7) = (-l) m y ^i-^,a) = \j -^Y lm {p,a). 

By exploiting (13), it is readily seen that isotropy (i.e., rotational invariance 
in law) entails 

oo I oo I 

Yl J2 ai m Yi m (x)=^2 Yl aimXlmigx) 

1=0 m=— I 1=0 m=— I 
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oo I I 

= Y Y a lm Y Dl m'm(9)Ylm>(x) 
1=0 m=— I m'=—l 



oo I / I \ 

Y Y [Y a imD l m , m {g)\Yi m ,(x), 

l=0m'=-l\ m =-l ) 



that is, 

(15) ( ai .)l(D l (g) ai .), 1 = 1,2,..., 

where the identity in law holds jointly over I and (aj.) denotes the 21 + 1 
vector (a;_ m , . . . , a; m ). We now recall the expression for the so-called Gaunt 
integral 



*l\m\ \-E) 1121712 \%) 1131713 {,•£) ">X 

(16) JS2 , 



= / (2Zi + l)(2Z 2 + l)(2Z 3 + l) (h h h\(h h h 
V 47T \0 J \mi 1712 ms 

which again requires a group-theoretical proof; indeed, the so-called Wigner's 
3j coefficients can be viewed as elements of the unitary matrices which al- 
ternate representation of the group of rotations; see Varshalovich, Moskalev 
and Khersonskii (1988), Vilenkin and Klymik (1991), Marinucci (2006) and 
Marinucci and Peccati (2007) for a more detailed discussion and explicit 
expressions. 

From (13) we have that under an arbitrary rotation the spherical har- 
monics transform as 

(17) ai m = Y D 7n>m(9)aim', 

in' 

and 

Ea hrni a hrn2 ai 3rn3 = Y ^mim'^^m^^mam^^Wi^m^m!, 

= Y D h ,{g)D h ,(g)D h , (a) 
' * minij"' 1712m 2 ya ' 7713177 .g va ' 

h h h\ ( h h h 



m'mLm., 



1 ) kihh 
h h h\ ( h h h 



/ \ m 1 m 2 rrii 



/ \m\ rri2 mz 

where we have used the identity 

h h h \ n h / \ n i 2 



bhhh-, 



^^TO 2 r7^g 



Y L1 , i2 f L3 f D h , (g)D l > , {g)D h , (g) 
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(18) 

\17li 771-2 rn 3. 

For a proof of (18), it suffices to use (16) and note that 



E 



u 



u 



h 



m\ m'n 



3 , )D h Ag)D h Ag)D h , (g) 

Tun I "iiwij " ' m^mj"' 77-13 m g V a > 



h h h 





(2Ji + l)(2Z 2 + l)(2Z 3 + l) 



f V D h ,(g)D l2 Ag)D h , (g)Y, _# Y, m , Y hm > 

Jo L / / / 



I'jm^mj 



r/;r 



h h h 




(2Zi + 1)(2Z 2 + 1)(2Z 3 + 1) 

47T 



S2 



[•Ml mi •* l2m,2 * 131713] &•£ 



= (h h h 

\m\ m2 1713 

For the first step to be justified, we need to ensure the Wigner's 3j coeffi- 
cients within the curly brackets is indeed different from zero; this condition 
is fulfilled provided h + h + h = even and the so-called triangle conditions 
are met, namely, li + lj < Ik for all choices of i,j, k = 1, 2, 3. 

We now turn to the issue of sample estimators. We can indeed show 
immediately that (4) is invariant to rotations; we have, in fact, 



1 



21 + 



I 53 \ a ^ 



1 



m=—l 



21 + 
1 



iS 



m=—l 



E^mm'W^ 



21 + 



j- 53 ai mi ai m2 53 Dl mmi (9)D l mm2 {g) 



1711,1712 



m=—l 



1 1 

"^ mi,iri2 m=-l 



2 
m| j 



in view of the orthonormality property Varshalovich, Moskalev and Khersonskii 

(1988) 

l 

/ j L'mmSd) "mm 2 \9) = "mi "h • 



m= — l 



A similar argument exploiting (18) and (17) shows indeed that the sample 
bispectrum is itself invariant to rotations. We refer again to Hu (2001), 
Marinucci (2006) and Marinucci (2008) for a much wider discussion and 
more properties. 
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A. 2. Some background on differential geometry. 

A. 2.1. Scalars, vectors and tensors. Our purpose is to establish intrinsic 
measures of the local curvature of a random field. Here, by intrinsic we 
mean "independent from the choice of coordinates," much the same way 
as our bispectrum statistics in the previous section. To clarify this issue, 
we recall some basic definition from differential geometry on the sphere; 
there exist many beautiful books on this issue [Adler and Taylor (2007), 
Bishop and Goldberg (1980)], and we refer to these for a proper account: we 
just report some basic facts to make the local curvature approach intuitive. 
Let M be a smooth manifold and write T x for the tangent plane at x £ M; 
we label T* the cotangent plane, that is, the dual space of T x [refer again 
to Adler and Taylor (2007) and Bishop and Goldberg (1980) for details and 
definitions] . Let (j) : M — ► K be some smooth function on the manifold; we 
say (j) is a, scalar function if it transforms under a change of coordinate g as 

4>{gx) = 4>(x) Vx £ M, that is, <fi := <ft o h~ . 

A rank one covariant tensor is a linear operator on the vector space T x with 
components {/»(■)} which transform according to the rule 

where we wrote y := {y 1 , . . . ,y n } = {gi(x), . . . ,g n (x)}. Likewise, a rank one 
contravariant tensor of dimension n is a linear operator {/i(-)}i=i,..., n whose 
components transform according to the rule 

i=l 

where we wrote as before y = gx. 

It is then possible to extend this definition to higher orders — for our pur- 
poses, rank 2 suffices. A rank two covariant tensor of dimension n is a 
bilinear operator {T uv } UtV= i t ^^ n whose components transform according to 
the rule 

(!9) T uv (y)=±T pq (x)^^. 

p,q=l " " 

Likewise, we can define contravariant and mixed rank two vectors, denoted 
by T uv and T" respectively. It is immediately seen from (19) that the usual 
Hessian matrix of a scalar function is a rank two covariant tensor. 

The previous concepts assume a nontrivial meaning in the presence of 
manifolds with a nonzero intrinsic curvature, such as the sphere. In such 
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spaces, we can introduce the metric tensor {g a b(-)} by imposing that the 
length of a vector X := {x 1 , . . . , x n } be given by 

n 
11*11? ■= E 9ab(x)x a X b . 
a,b=l 

Thus, {g a b( x )}a,b=i,...,n is a rank two tensor. Some examples: 

• Euclidean case: we have 

n 

\\X\\e ■= J2 ^^a ^at is, g ab = 5 b a . 

a,b=l 

• The two sphere S 2 : in spherical coordinates, for f2 := {$,99} the usual 
metric is given by 

(20) {gab(n)} a , b= ^=(l J 2 ^j. 

The contravariant or reciprocal metric tensor {g {x)} a: b=i,..., n is defined 
by the requirement that 

Y,9 ab (x)9bc(x) = S a c , 

b 

that is, it represents the elements of the inverse matrix g . By means of 
this tensor operator, we can define more generally contravariant vectors, 
denoting, for instance, 

T a :=Y,9abT b and consequently, T a = £>"% 
b b 

so that we ensure invariance, that is, 

||{T4 a=1 ,...,j| 5 = ^r a r a = ^ 5a6 r a T b = ^^T a r 6 =:||{r a }a=i,...,n|| 5 . 

1 a,b a,b 

Covariant tensors can be likewise introduced, that is, 

rpa \ ■* „acrp rpab \ * „acrpb \ * ^ac Jbdrp 

I b=2-^9 -Lac, ± = 2^9 I c=2-^9 9 +cd- 

c c c ,d 

Let us now investigate the behavior of second order derivatives under coor- 
dinate transformations. We have 

d 2 (f) d ( d(j) 



QyU QyV QyU ^ QyV 



J^ dx p dx q J^ d 2 x p 

p,q=l a a p=l a w 
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It is therefore clear that for nonlinear coordinate transformations (i.e., such 
that {d 2 x p /dy u dy v } 7^ 0) standard second order derivatives do not act as 
a rank two tensor, that is, they depend on the coordinate choice and do 
not represent intrinsic quantities. To overcome this problem, we need to 
introduce the well-known Christoffel coefficients T ik , which, by assumption, 
satisfy the transformation laws 

— i y^ r n ® x ^ d 00 ™ 1 ®y l v^ d 2 ^ dy l 

ik = 2^ jm Q~i q k dxn + 2^ q i q k Q x j ■ 
m,n,j a a 3 

It is then easy to check that 

I rn,n,j a a j a a 

and, hence, 

dx^ dx r 



J i,k 



2~2 T ik( L h<i>)l = J2 ^< m ~ J2 Tl jm<Pn 1 
I m,j \ n J 



Qyl Qyk ' 



that is, {(pj, m — Sn^jm^n} is actually a rank two tensor. The previous 
discussion motivates the following: 

Definition. The covariant derivative of the rank one tensor Tj is given 
by 

Ti-k '-T^k -^2^ik T h 
i 

where Tj^ = dTi/dx denotes standard derivative. 

Let us now specify the previous definition to the sphere. In terms of the 
metric tensor, the Christoffel symbols can be shown to be equal to 



( 21 ) T \j = 2~L\Y^ki,j + 9kj,i ~ 9ij,k)>- 



For instance, in the Euclidean case gkij = gkj,i = 9ij,k = 0, hence, T\ ■ = 0; 
thus, covariant derivatives coincide with standard calculus operators. On the 
contrary, for the sphere S 2 , we have g#@ = g = 1, g w = sin 2 -d = (g w ) _1 , 
9-dtp = 9 = 0. Hence, we obtain 

i? = 5i?i?,v5 = 9&<p,$ = 9$<p,<p = 9<p<p,<p = 
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and gipu,^ = 2 sin?? cos??; to summarize, the Christoffel symbols on the sphere 
are 

1 W — L &0 — u i 

gfP I 

T % = if{gw,tf + g^-g^}= ■ 2 2 sin?? cos?? = cot??, 
^ _■ sin 1/ 

m &<p 

^tp-a = —^-{giptifl + am^ - gtp&,&} + —^-{g vv ,-a + g^ v , v - g<p#,<p} = 0, 

„w 

1 (pep ~ ~^~g<pip,tp — u i 

g m 

r w = -o~{9tp0,tp + g<p$,<p - g<pip,&} = - sin?? cos??. 

Hence, we obtain the following results for the covariant derivatives of the 
spherical harmonics: 

llm;-d(p — ilmfiip J- </,#-* fan,</S J- tp-Q^lm,^ = ±lm,flip c Crt ""lm,(pi 
Xlm;tptp = Ilm,(ptp -I tpip^lm,<p ^ <^j-Mm,$ — *lm,(p(p t Sin?/ COS?/i lm,-&- 

The previous expressions provide the clue for the computation of the bilinear 
form 

(22") H* ■= C^' m T;tiip\ _ {J2lm a lmYl m ;m J2lm a lmXlm;-d(p\ 

\-L;$tp -L;tpip/ \ 2-^lm ®lm *lm;'&(p L^tlm ®lm *lm;tpip J 

To obtain (12), we need to introduce a final, quite subtle point. (22) defines 
a bilinear form H* : (T x T) -> 1 acting on the tensor product of the tangent 
plane with itself; in order to be able to evaluate consistently the eigenvalues, 
we must transform this form into the corresponding linear operator H : T — ► 
T, where H := g~ l H* [actually we considered the symmetrized expression 
H := g ' H*g ' , where g denotes as before the metric tensor on the 
sphere, see (20)]. This explains the appearance of the sin?? factors at the 
denominators in (12) — we refer again to Bishop and Goldberg (1980) for 
more details and explanations. 
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